// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SUPERLUSUPPORT_H
#define EIGEN_SUPERLUSUPPORT_H

namespace Eigen {

#if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
#define DECL_GSSVX(PREFIX, FLOATTYPE, KEYTYPE)                                                                         \
	extern "C"                                                                                                         \
	{                                                                                                                  \
		extern void PREFIX##gssvx(superlu_options_t*,                                                                  \
								  SuperMatrix*,                                                                        \
								  int*,                                                                                \
								  int*,                                                                                \
								  int*,                                                                                \
								  char*,                                                                               \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  void*,                                                                               \
								  int,                                                                                 \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  GlobalLU_t*,                                                                         \
								  mem_usage_t*,                                                                        \
								  SuperLUStat_t*,                                                                      \
								  int*);                                                                               \
	}                                                                                                                  \
	inline float SuperLU_gssvx(superlu_options_t* options,                                                             \
							   SuperMatrix* A,                                                                         \
							   int* perm_c,                                                                            \
							   int* perm_r,                                                                            \
							   int* etree,                                                                             \
							   char* equed,                                                                            \
							   FLOATTYPE* R,                                                                           \
							   FLOATTYPE* C,                                                                           \
							   SuperMatrix* L,                                                                         \
							   SuperMatrix* U,                                                                         \
							   void* work,                                                                             \
							   int lwork,                                                                              \
							   SuperMatrix* B,                                                                         \
							   SuperMatrix* X,                                                                         \
							   FLOATTYPE* recip_pivot_growth,                                                          \
							   FLOATTYPE* rcond,                                                                       \
							   FLOATTYPE* ferr,                                                                        \
							   FLOATTYPE* berr,                                                                        \
							   SuperLUStat_t* stats,                                                                   \
							   int* info,                                                                              \
							   KEYTYPE)                                                                                \
	{                                                                                                                  \
		mem_usage_t mem_usage;                                                                                         \
		GlobalLU_t gLU;                                                                                                \
		PREFIX##gssvx(options,                                                                                         \
					  A,                                                                                               \
					  perm_c,                                                                                          \
					  perm_r,                                                                                          \
					  etree,                                                                                           \
					  equed,                                                                                           \
					  R,                                                                                               \
					  C,                                                                                               \
					  L,                                                                                               \
					  U,                                                                                               \
					  work,                                                                                            \
					  lwork,                                                                                           \
					  B,                                                                                               \
					  X,                                                                                               \
					  recip_pivot_growth,                                                                              \
					  rcond,                                                                                           \
					  ferr,                                                                                            \
					  berr,                                                                                            \
					  &gLU,                                                                                            \
					  &mem_usage,                                                                                      \
					  stats,                                                                                           \
					  info);                                                                                           \
		return mem_usage.for_lu; /* bytes used by the factor storage */                                                \
	}
#else // version < 5.0
#define DECL_GSSVX(PREFIX, FLOATTYPE, KEYTYPE)                                                                         \
	extern "C"                                                                                                         \
	{                                                                                                                  \
		extern void PREFIX##gssvx(superlu_options_t*,                                                                  \
								  SuperMatrix*,                                                                        \
								  int*,                                                                                \
								  int*,                                                                                \
								  int*,                                                                                \
								  char*,                                                                               \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  void*,                                                                               \
								  int,                                                                                 \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  mem_usage_t*,                                                                        \
								  SuperLUStat_t*,                                                                      \
								  int*);                                                                               \
	}                                                                                                                  \
	inline float SuperLU_gssvx(superlu_options_t* options,                                                             \
							   SuperMatrix* A,                                                                         \
							   int* perm_c,                                                                            \
							   int* perm_r,                                                                            \
							   int* etree,                                                                             \
							   char* equed,                                                                            \
							   FLOATTYPE* R,                                                                           \
							   FLOATTYPE* C,                                                                           \
							   SuperMatrix* L,                                                                         \
							   SuperMatrix* U,                                                                         \
							   void* work,                                                                             \
							   int lwork,                                                                              \
							   SuperMatrix* B,                                                                         \
							   SuperMatrix* X,                                                                         \
							   FLOATTYPE* recip_pivot_growth,                                                          \
							   FLOATTYPE* rcond,                                                                       \
							   FLOATTYPE* ferr,                                                                        \
							   FLOATTYPE* berr,                                                                        \
							   SuperLUStat_t* stats,                                                                   \
							   int* info,                                                                              \
							   KEYTYPE)                                                                                \
	{                                                                                                                  \
		mem_usage_t mem_usage;                                                                                         \
		PREFIX##gssvx(options,                                                                                         \
					  A,                                                                                               \
					  perm_c,                                                                                          \
					  perm_r,                                                                                          \
					  etree,                                                                                           \
					  equed,                                                                                           \
					  R,                                                                                               \
					  C,                                                                                               \
					  L,                                                                                               \
					  U,                                                                                               \
					  work,                                                                                            \
					  lwork,                                                                                           \
					  B,                                                                                               \
					  X,                                                                                               \
					  recip_pivot_growth,                                                                              \
					  rcond,                                                                                           \
					  ferr,                                                                                            \
					  berr,                                                                                            \
					  &mem_usage,                                                                                      \
					  stats,                                                                                           \
					  info);                                                                                           \
		return mem_usage.for_lu; /* bytes used by the factor storage */                                                \
	}
#endif

DECL_GSSVX(s, float, float)
DECL_GSSVX(c, float, std::complex<float>)
DECL_GSSVX(d, double, double)
DECL_GSSVX(z, double, std::complex<double>)

#ifdef MILU_ALPHA
#define EIGEN_SUPERLU_HAS_ILU
#endif

#ifdef EIGEN_SUPERLU_HAS_ILU

// similarly for the incomplete factorization using gsisx
#define DECL_GSISX(PREFIX, FLOATTYPE, KEYTYPE)                                                                         \
	extern "C"                                                                                                         \
	{                                                                                                                  \
		extern void PREFIX##gsisx(superlu_options_t*,                                                                  \
								  SuperMatrix*,                                                                        \
								  int*,                                                                                \
								  int*,                                                                                \
								  int*,                                                                                \
								  char*,                                                                               \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  void*,                                                                               \
								  int,                                                                                 \
								  SuperMatrix*,                                                                        \
								  SuperMatrix*,                                                                        \
								  FLOATTYPE*,                                                                          \
								  FLOATTYPE*,                                                                          \
								  mem_usage_t*,                                                                        \
								  SuperLUStat_t*,                                                                      \
								  int*);                                                                               \
	}                                                                                                                  \
	inline float SuperLU_gsisx(superlu_options_t* options,                                                             \
							   SuperMatrix* A,                                                                         \
							   int* perm_c,                                                                            \
							   int* perm_r,                                                                            \
							   int* etree,                                                                             \
							   char* equed,                                                                            \
							   FLOATTYPE* R,                                                                           \
							   FLOATTYPE* C,                                                                           \
							   SuperMatrix* L,                                                                         \
							   SuperMatrix* U,                                                                         \
							   void* work,                                                                             \
							   int lwork,                                                                              \
							   SuperMatrix* B,                                                                         \
							   SuperMatrix* X,                                                                         \
							   FLOATTYPE* recip_pivot_growth,                                                          \
							   FLOATTYPE* rcond,                                                                       \
							   SuperLUStat_t* stats,                                                                   \
							   int* info,                                                                              \
							   KEYTYPE)                                                                                \
	{                                                                                                                  \
		mem_usage_t mem_usage;                                                                                         \
		PREFIX##gsisx(options,                                                                                         \
					  A,                                                                                               \
					  perm_c,                                                                                          \
					  perm_r,                                                                                          \
					  etree,                                                                                           \
					  equed,                                                                                           \
					  R,                                                                                               \
					  C,                                                                                               \
					  L,                                                                                               \
					  U,                                                                                               \
					  work,                                                                                            \
					  lwork,                                                                                           \
					  B,                                                                                               \
					  X,                                                                                               \
					  recip_pivot_growth,                                                                              \
					  rcond,                                                                                           \
					  &mem_usage,                                                                                      \
					  stats,                                                                                           \
					  info);                                                                                           \
		return mem_usage.for_lu; /* bytes used by the factor storage */                                                \
	}

DECL_GSISX(s, float, float)
DECL_GSISX(c, float, std::complex<float>)
DECL_GSISX(d, double, double)
DECL_GSISX(z, double, std::complex<double>)

#endif

template<typename MatrixType>
struct SluMatrixMapHelper;

/** \internal
 *
 * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
 * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
 *
 * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
 */
struct SluMatrix : SuperMatrix
{
	SluMatrix() { Store = &storage; }

	SluMatrix(const SluMatrix& other)
		: SuperMatrix(other)
	{
		Store = &storage;
		storage = other.storage;
	}

	SluMatrix& operator=(const SluMatrix& other)
	{
		SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
		Store = &storage;
		storage = other.storage;
		return *this;
	}

	struct
	{
		union
		{
			int nnz;
			int lda;
		};
		void* values;
		int* innerInd;
		int* outerInd;
	} storage;

	void setStorageType(Stype_t t)
	{
		Stype = t;
		if (t == SLU_NC || t == SLU_NR || t == SLU_DN)
			Store = &storage;
		else {
			eigen_assert(false && "storage type not supported");
			Store = 0;
		}
	}

	template<typename Scalar>
	void setScalarType()
	{
		if (internal::is_same<Scalar, float>::value)
			Dtype = SLU_S;
		else if (internal::is_same<Scalar, double>::value)
			Dtype = SLU_D;
		else if (internal::is_same<Scalar, std::complex<float>>::value)
			Dtype = SLU_C;
		else if (internal::is_same<Scalar, std::complex<double>>::value)
			Dtype = SLU_Z;
		else {
			eigen_assert(false && "Scalar type not supported by SuperLU");
		}
	}

	template<typename MatrixType>
	static SluMatrix Map(MatrixBase<MatrixType>& _mat)
	{
		MatrixType& mat(_mat.derived());
		eigen_assert(((MatrixType::Flags & RowMajorBit) != RowMajorBit) &&
					 "row-major dense matrices are not supported by SuperLU");
		SluMatrix res;
		res.setStorageType(SLU_DN);
		res.setScalarType<typename MatrixType::Scalar>();
		res.Mtype = SLU_GE;

		res.nrow = internal::convert_index<int>(mat.rows());
		res.ncol = internal::convert_index<int>(mat.cols());

		res.storage.lda =
			internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
		res.storage.values = (void*)(mat.data());
		return res;
	}

	template<typename MatrixType>
	static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
	{
		MatrixType& mat(a_mat.derived());
		SluMatrix res;
		if ((MatrixType::Flags & RowMajorBit) == RowMajorBit) {
			res.setStorageType(SLU_NR);
			res.nrow = internal::convert_index<int>(mat.cols());
			res.ncol = internal::convert_index<int>(mat.rows());
		} else {
			res.setStorageType(SLU_NC);
			res.nrow = internal::convert_index<int>(mat.rows());
			res.ncol = internal::convert_index<int>(mat.cols());
		}

		res.Mtype = SLU_GE;

		res.storage.nnz = internal::convert_index<int>(mat.nonZeros());
		res.storage.values = mat.valuePtr();
		res.storage.innerInd = mat.innerIndexPtr();
		res.storage.outerInd = mat.outerIndexPtr();

		res.setScalarType<typename MatrixType::Scalar>();

		// FIXME the following is not very accurate
		if (int(MatrixType::Flags) & int(Upper))
			res.Mtype = SLU_TRU;
		if (int(MatrixType::Flags) & int(Lower))
			res.Mtype = SLU_TRL;

		eigen_assert(((int(MatrixType::Flags) & int(SelfAdjoint)) == 0) &&
					 "SelfAdjoint matrix shape not supported by SuperLU");

		return res;
	}
};

template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
struct SluMatrixMapHelper<Matrix<Scalar, Rows, Cols, Options, MRows, MCols>>
{
	typedef Matrix<Scalar, Rows, Cols, Options, MRows, MCols> MatrixType;
	static void run(MatrixType& mat, SluMatrix& res)
	{
		eigen_assert(((Options & RowMajor) != RowMajor) && "row-major dense matrices is not supported by SuperLU");
		res.setStorageType(SLU_DN);
		res.setScalarType<Scalar>();
		res.Mtype = SLU_GE;

		res.nrow = mat.rows();
		res.ncol = mat.cols();

		res.storage.lda = mat.outerStride();
		res.storage.values = mat.data();
	}
};

template<typename Derived>
struct SluMatrixMapHelper<SparseMatrixBase<Derived>>
{
	typedef Derived MatrixType;
	static void run(MatrixType& mat, SluMatrix& res)
	{
		if ((MatrixType::Flags & RowMajorBit) == RowMajorBit) {
			res.setStorageType(SLU_NR);
			res.nrow = mat.cols();
			res.ncol = mat.rows();
		} else {
			res.setStorageType(SLU_NC);
			res.nrow = mat.rows();
			res.ncol = mat.cols();
		}

		res.Mtype = SLU_GE;

		res.storage.nnz = mat.nonZeros();
		res.storage.values = mat.valuePtr();
		res.storage.innerInd = mat.innerIndexPtr();
		res.storage.outerInd = mat.outerIndexPtr();

		res.setScalarType<typename MatrixType::Scalar>();

		// FIXME the following is not very accurate
		if (MatrixType::Flags & Upper)
			res.Mtype = SLU_TRU;
		if (MatrixType::Flags & Lower)
			res.Mtype = SLU_TRL;

		eigen_assert(((MatrixType::Flags & SelfAdjoint) == 0) && "SelfAdjoint matrix shape not supported by SuperLU");
	}
};

namespace internal {

template<typename MatrixType>
SluMatrix
asSluMatrix(MatrixType& mat)
{
	return SluMatrix::Map(mat);
}

/** View a Super LU matrix as an Eigen expression */
template<typename Scalar, int Flags, typename Index>
MappedSparseMatrix<Scalar, Flags, Index>
map_superlu(SluMatrix& sluMat)
{
	eigen_assert(((Flags & RowMajor) == RowMajor && sluMat.Stype == SLU_NR) ||
				 ((Flags & ColMajor) == ColMajor && sluMat.Stype == SLU_NC));

	Index outerSize = (Flags & RowMajor) == RowMajor ? sluMat.ncol : sluMat.nrow;

	return MappedSparseMatrix<Scalar, Flags, Index>(sluMat.nrow,
													sluMat.ncol,
													sluMat.storage.outerInd[outerSize],
													sluMat.storage.outerInd,
													sluMat.storage.innerInd,
													reinterpret_cast<Scalar*>(sluMat.storage.values));
}

} // end namespace internal

/** \ingroup SuperLUSupport_Module
 * \class SuperLUBase
 * \brief The base class for the direct and incomplete LU factorization of SuperLU
 */
template<typename _MatrixType, typename Derived>
class SuperLUBase : public SparseSolverBase<Derived>
{
  protected:
	typedef SparseSolverBase<Derived> Base;
	using Base::derived;
	using Base::m_isInitialized;

  public:
	typedef _MatrixType MatrixType;
	typedef typename MatrixType::Scalar Scalar;
	typedef typename MatrixType::RealScalar RealScalar;
	typedef typename MatrixType::StorageIndex StorageIndex;
	typedef Matrix<Scalar, Dynamic, 1> Vector;
	typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
	typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
	typedef Map<PermutationMatrix<Dynamic, Dynamic, int>> PermutationMap;
	typedef SparseMatrix<Scalar> LUMatrixType;
	enum
	{
		ColsAtCompileTime = MatrixType::ColsAtCompileTime,
		MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
	};

  public:
	SuperLUBase() {}

	~SuperLUBase() { clearFactors(); }

	inline Index rows() const { return m_matrix.rows(); }
	inline Index cols() const { return m_matrix.cols(); }

	/** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
	inline superlu_options_t& options() { return m_sluOptions; }

	/** \brief Reports whether previous computation was successful.
	 *
	 * \returns \c Success if computation was successful,
	 *          \c NumericalIssue if the matrix.appears to be negative.
	 */
	ComputationInfo info() const
	{
		eigen_assert(m_isInitialized && "Decomposition is not initialized.");
		return m_info;
	}

	/** Computes the sparse Cholesky decomposition of \a matrix */
	void compute(const MatrixType& matrix)
	{
		derived().analyzePattern(matrix);
		derived().factorize(matrix);
	}

	/** Performs a symbolic decomposition on the sparcity of \a matrix.
	 *
	 * This function is particularly useful when solving for several problems having the same structure.
	 *
	 * \sa factorize()
	 */
	void analyzePattern(const MatrixType& /*matrix*/)
	{
		m_isInitialized = true;
		m_info = Success;
		m_analysisIsOk = true;
		m_factorizationIsOk = false;
	}

	template<typename Stream>
	void dumpMemory(Stream& /*s*/)
	{
	}

  protected:
	void initFactorization(const MatrixType& a)
	{
		set_default_options(&this->m_sluOptions);

		const Index size = a.rows();
		m_matrix = a;

		m_sluA = internal::asSluMatrix(m_matrix);
		clearFactors();

		m_p.resize(size);
		m_q.resize(size);
		m_sluRscale.resize(size);
		m_sluCscale.resize(size);
		m_sluEtree.resize(size);

		// set empty B and X
		m_sluB.setStorageType(SLU_DN);
		m_sluB.setScalarType<Scalar>();
		m_sluB.Mtype = SLU_GE;
		m_sluB.storage.values = 0;
		m_sluB.nrow = 0;
		m_sluB.ncol = 0;
		m_sluB.storage.lda = internal::convert_index<int>(size);
		m_sluX = m_sluB;

		m_extractedDataAreDirty = true;
	}

	void init()
	{
		m_info = InvalidInput;
		m_isInitialized = false;
		m_sluL.Store = 0;
		m_sluU.Store = 0;
	}

	void extractData() const;

	void clearFactors()
	{
		if (m_sluL.Store)
			Destroy_SuperNode_Matrix(&m_sluL);
		if (m_sluU.Store)
			Destroy_CompCol_Matrix(&m_sluU);

		m_sluL.Store = 0;
		m_sluU.Store = 0;

		memset(&m_sluL, 0, sizeof m_sluL);
		memset(&m_sluU, 0, sizeof m_sluU);
	}

	// cached data to reduce reallocation, etc.
	mutable LUMatrixType m_l;
	mutable LUMatrixType m_u;
	mutable IntColVectorType m_p;
	mutable IntRowVectorType m_q;

	mutable LUMatrixType m_matrix; // copy of the factorized matrix
	mutable SluMatrix m_sluA;
	mutable SuperMatrix m_sluL, m_sluU;
	mutable SluMatrix m_sluB, m_sluX;
	mutable SuperLUStat_t m_sluStat;
	mutable superlu_options_t m_sluOptions;
	mutable std::vector<int> m_sluEtree;
	mutable Matrix<RealScalar, Dynamic, 1> m_sluRscale, m_sluCscale;
	mutable Matrix<RealScalar, Dynamic, 1> m_sluFerr, m_sluBerr;
	mutable char m_sluEqued;

	mutable ComputationInfo m_info;
	int m_factorizationIsOk;
	int m_analysisIsOk;
	mutable bool m_extractedDataAreDirty;

  private:
	SuperLUBase(SuperLUBase&) {}
};

/** \ingroup SuperLUSupport_Module
 * \class SuperLU
 * \brief A sparse direct LU factorization and solver based on the SuperLU library
 *
 * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
 * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
 * X and B can be either dense or sparse.
 *
 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 *
 * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
 *
 * \implsparsesolverconcept
 *
 * \sa \ref TutorialSparseSolverConcept, class SparseLU
 */
template<typename _MatrixType>
class SuperLU : public SuperLUBase<_MatrixType, SuperLU<_MatrixType>>
{
  public:
	typedef SuperLUBase<_MatrixType, SuperLU> Base;
	typedef _MatrixType MatrixType;
	typedef typename Base::Scalar Scalar;
	typedef typename Base::RealScalar RealScalar;
	typedef typename Base::StorageIndex StorageIndex;
	typedef typename Base::IntRowVectorType IntRowVectorType;
	typedef typename Base::IntColVectorType IntColVectorType;
	typedef typename Base::PermutationMap PermutationMap;
	typedef typename Base::LUMatrixType LUMatrixType;
	typedef TriangularView<LUMatrixType, Lower | UnitDiag> LMatrixType;
	typedef TriangularView<LUMatrixType, Upper> UMatrixType;

  public:
	using Base::_solve_impl;

	SuperLU()
		: Base()
	{
		init();
	}

	explicit SuperLU(const MatrixType& matrix)
		: Base()
	{
		init();
		Base::compute(matrix);
	}

	~SuperLU() {}

	/** Performs a symbolic decomposition on the sparcity of \a matrix.
	 *
	 * This function is particularly useful when solving for several problems having the same structure.
	 *
	 * \sa factorize()
	 */
	void analyzePattern(const MatrixType& matrix)
	{
		m_info = InvalidInput;
		m_isInitialized = false;
		Base::analyzePattern(matrix);
	}

	/** Performs a numeric decomposition of \a matrix
	 *
	 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been
	 * performed.
	 *
	 * \sa analyzePattern()
	 */
	void factorize(const MatrixType& matrix);

	/** \internal */
	template<typename Rhs, typename Dest>
	void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const;

	inline const LMatrixType& matrixL() const
	{
		if (m_extractedDataAreDirty)
			this->extractData();
		return m_l;
	}

	inline const UMatrixType& matrixU() const
	{
		if (m_extractedDataAreDirty)
			this->extractData();
		return m_u;
	}

	inline const IntColVectorType& permutationP() const
	{
		if (m_extractedDataAreDirty)
			this->extractData();
		return m_p;
	}

	inline const IntRowVectorType& permutationQ() const
	{
		if (m_extractedDataAreDirty)
			this->extractData();
		return m_q;
	}

	Scalar determinant() const;

  protected:
	using Base::m_l;
	using Base::m_matrix;
	using Base::m_p;
	using Base::m_q;
	using Base::m_sluA;
	using Base::m_sluB;
	using Base::m_sluBerr;
	using Base::m_sluCscale;
	using Base::m_sluEqued;
	using Base::m_sluEtree;
	using Base::m_sluFerr;
	using Base::m_sluL;
	using Base::m_sluOptions;
	using Base::m_sluRscale;
	using Base::m_sluStat;
	using Base::m_sluU;
	using Base::m_sluX;
	using Base::m_u;

	using Base::m_analysisIsOk;
	using Base::m_extractedDataAreDirty;
	using Base::m_factorizationIsOk;
	using Base::m_info;
	using Base::m_isInitialized;

	void init()
	{
		Base::init();

		set_default_options(&this->m_sluOptions);
		m_sluOptions.PrintStat = NO;
		m_sluOptions.ConditionNumber = NO;
		m_sluOptions.Trans = NOTRANS;
		m_sluOptions.ColPerm = COLAMD;
	}

  private:
	SuperLU(SuperLU&) {}
};

template<typename MatrixType>
void
SuperLU<MatrixType>::factorize(const MatrixType& a)
{
	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
	if (!m_analysisIsOk) {
		m_info = InvalidInput;
		return;
	}

	this->initFactorization(a);

	m_sluOptions.ColPerm = COLAMD;
	int info = 0;
	RealScalar recip_pivot_growth, rcond;
	RealScalar ferr, berr;

	StatInit(&m_sluStat);
	SuperLU_gssvx(&m_sluOptions,
				  &m_sluA,
				  m_q.data(),
				  m_p.data(),
				  &m_sluEtree[0],
				  &m_sluEqued,
				  &m_sluRscale[0],
				  &m_sluCscale[0],
				  &m_sluL,
				  &m_sluU,
				  NULL,
				  0,
				  &m_sluB,
				  &m_sluX,
				  &recip_pivot_growth,
				  &rcond,
				  &ferr,
				  &berr,
				  &m_sluStat,
				  &info,
				  Scalar());
	StatFree(&m_sluStat);

	m_extractedDataAreDirty = true;

	// FIXME how to better check for errors ???
	m_info = info == 0 ? Success : NumericalIssue;
	m_factorizationIsOk = true;
}

template<typename MatrixType>
template<typename Rhs, typename Dest>
void
SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& x) const
{
	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call "
										"either compute() or analyzePattern()/factorize()");

	const Index rhsCols = b.cols();
	eigen_assert(m_matrix.rows() == b.rows());

	m_sluOptions.Trans = NOTRANS;
	m_sluOptions.Fact = FACTORED;
	m_sluOptions.IterRefine = NOREFINE;

	m_sluFerr.resize(rhsCols);
	m_sluBerr.resize(rhsCols);

	Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor>> b_ref(b);
	Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor>> x_ref(x);

	m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
	m_sluX = SluMatrix::Map(x_ref.const_cast_derived());

	typename Rhs::PlainObject b_cpy;
	if (m_sluEqued != 'N') {
		b_cpy = b;
		m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
	}

	StatInit(&m_sluStat);
	int info = 0;
	RealScalar recip_pivot_growth, rcond;
	SuperLU_gssvx(&m_sluOptions,
				  &m_sluA,
				  m_q.data(),
				  m_p.data(),
				  &m_sluEtree[0],
				  &m_sluEqued,
				  &m_sluRscale[0],
				  &m_sluCscale[0],
				  &m_sluL,
				  &m_sluU,
				  NULL,
				  0,
				  &m_sluB,
				  &m_sluX,
				  &recip_pivot_growth,
				  &rcond,
				  &m_sluFerr[0],
				  &m_sluBerr[0],
				  &m_sluStat,
				  &info,
				  Scalar());
	StatFree(&m_sluStat);

	if (x.derived().data() != x_ref.data())
		x = x_ref;

	m_info = info == 0 ? Success : NumericalIssue;
}

// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
//
//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
//
//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
//
template<typename MatrixType, typename Derived>
void
SuperLUBase<MatrixType, Derived>::extractData() const
{
	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must "
										"first call either compute() or analyzePattern()/factorize()");
	if (m_extractedDataAreDirty) {
		int upper;
		int fsupc, istart, nsupr;
		int lastl = 0, lastu = 0;
		SCformat* Lstore = static_cast<SCformat*>(m_sluL.Store);
		NCformat* Ustore = static_cast<NCformat*>(m_sluU.Store);
		Scalar* SNptr;

		const Index size = m_matrix.rows();
		m_l.resize(size, size);
		m_l.resizeNonZeros(Lstore->nnz);
		m_u.resize(size, size);
		m_u.resizeNonZeros(Ustore->nnz);

		int* Lcol = m_l.outerIndexPtr();
		int* Lrow = m_l.innerIndexPtr();
		Scalar* Lval = m_l.valuePtr();

		int* Ucol = m_u.outerIndexPtr();
		int* Urow = m_u.innerIndexPtr();
		Scalar* Uval = m_u.valuePtr();

		Ucol[0] = 0;
		Ucol[0] = 0;

		/* for each supernode */
		for (int k = 0; k <= Lstore->nsuper; ++k) {
			fsupc = L_FST_SUPC(k);
			istart = L_SUB_START(fsupc);
			nsupr = L_SUB_START(fsupc + 1) - istart;
			upper = 1;

			/* for each column in the supernode */
			for (int j = fsupc; j < L_FST_SUPC(k + 1); ++j) {
				SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];

				/* Extract U */
				for (int i = U_NZ_START(j); i < U_NZ_START(j + 1); ++i) {
					Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
					/* Matlab doesn't like explicit zero. */
					if (Uval[lastu] != 0.0)
						Urow[lastu++] = U_SUB(i);
				}
				for (int i = 0; i < upper; ++i) {
					/* upper triangle in the supernode */
					Uval[lastu] = SNptr[i];
					/* Matlab doesn't like explicit zero. */
					if (Uval[lastu] != 0.0)
						Urow[lastu++] = L_SUB(istart + i);
				}
				Ucol[j + 1] = lastu;

				/* Extract L */
				Lval[lastl] = 1.0; /* unit diagonal */
				Lrow[lastl++] = L_SUB(istart + upper - 1);
				for (int i = upper; i < nsupr; ++i) {
					Lval[lastl] = SNptr[i];
					/* Matlab doesn't like explicit zero. */
					if (Lval[lastl] != 0.0)
						Lrow[lastl++] = L_SUB(istart + i);
				}
				Lcol[j + 1] = lastl;

				++upper;
			} /* for j ... */

		} /* for k ... */

		// squeeze the matrices :
		m_l.resizeNonZeros(lastl);
		m_u.resizeNonZeros(lastu);

		m_extractedDataAreDirty = false;
	}
}

template<typename MatrixType>
typename SuperLU<MatrixType>::Scalar
SuperLU<MatrixType>::determinant() const
{
	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you "
										"must first call either compute() or analyzePattern()/factorize()");

	if (m_extractedDataAreDirty)
		this->extractData();

	Scalar det = Scalar(1);
	for (int j = 0; j < m_u.cols(); ++j) {
		if (m_u.outerIndexPtr()[j + 1] - m_u.outerIndexPtr()[j] > 0) {
			int lastId = m_u.outerIndexPtr()[j + 1] - 1;
			eigen_assert(m_u.innerIndexPtr()[lastId] <= j);
			if (m_u.innerIndexPtr()[lastId] == j)
				det *= m_u.valuePtr()[lastId];
		}
	}
	if (PermutationMap(m_p.data(), m_p.size()).determinant() * PermutationMap(m_q.data(), m_q.size()).determinant() < 0)
		det = -det;
	if (m_sluEqued != 'N')
		return det / m_sluRscale.prod() / m_sluCscale.prod();
	else
		return det;
}

#ifdef EIGEN_PARSED_BY_DOXYGEN
#define EIGEN_SUPERLU_HAS_ILU
#endif

#ifdef EIGEN_SUPERLU_HAS_ILU

/** \ingroup SuperLUSupport_Module
 * \class SuperILU
 * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
 *
 * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU
 * factorization using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear
 * solvers.
 *
 * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
 *
 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
 *
 * \implsparsesolverconcept
 *
 * \sa \ref TutorialSparseSolverConcept, class IncompleteLUT, class ConjugateGradient, class BiCGSTAB
 */

template<typename _MatrixType>
class SuperILU : public SuperLUBase<_MatrixType, SuperILU<_MatrixType>>
{
  public:
	typedef SuperLUBase<_MatrixType, SuperILU> Base;
	typedef _MatrixType MatrixType;
	typedef typename Base::Scalar Scalar;
	typedef typename Base::RealScalar RealScalar;

  public:
	using Base::_solve_impl;

	SuperILU()
		: Base()
	{
		init();
	}

	SuperILU(const MatrixType& matrix)
		: Base()
	{
		init();
		Base::compute(matrix);
	}

	~SuperILU() {}

	/** Performs a symbolic decomposition on the sparcity of \a matrix.
	 *
	 * This function is particularly useful when solving for several problems having the same structure.
	 *
	 * \sa factorize()
	 */
	void analyzePattern(const MatrixType& matrix) { Base::analyzePattern(matrix); }

	/** Performs a numeric decomposition of \a matrix
	 *
	 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been
	 * performed.
	 *
	 * \sa analyzePattern()
	 */
	void factorize(const MatrixType& matrix);

#ifndef EIGEN_PARSED_BY_DOXYGEN
	/** \internal */
	template<typename Rhs, typename Dest>
	void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const;
#endif // EIGEN_PARSED_BY_DOXYGEN

  protected:
	using Base::m_l;
	using Base::m_matrix;
	using Base::m_p;
	using Base::m_q;
	using Base::m_sluA;
	using Base::m_sluB;
	using Base::m_sluBerr;
	using Base::m_sluCscale;
	using Base::m_sluEqued;
	using Base::m_sluEtree;
	using Base::m_sluFerr;
	using Base::m_sluL;
	using Base::m_sluOptions;
	using Base::m_sluRscale;
	using Base::m_sluStat;
	using Base::m_sluU;
	using Base::m_sluX;
	using Base::m_u;

	using Base::m_analysisIsOk;
	using Base::m_extractedDataAreDirty;
	using Base::m_factorizationIsOk;
	using Base::m_info;
	using Base::m_isInitialized;

	void init()
	{
		Base::init();

		ilu_set_default_options(&m_sluOptions);
		m_sluOptions.PrintStat = NO;
		m_sluOptions.ConditionNumber = NO;
		m_sluOptions.Trans = NOTRANS;
		m_sluOptions.ColPerm = MMD_AT_PLUS_A;

		// no attempt to preserve column sum
		m_sluOptions.ILU_MILU = SILU;
		// only basic ILU(k) support -- no direct control over memory consumption
		// better to use ILU_DropRule = DROP_BASIC | DROP_AREA
		// and set ILU_FillFactor to max memory growth
		m_sluOptions.ILU_DropRule = DROP_BASIC;
		m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision() * 10;
	}

  private:
	SuperILU(SuperILU&) {}
};

template<typename MatrixType>
void
SuperILU<MatrixType>::factorize(const MatrixType& a)
{
	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
	if (!m_analysisIsOk) {
		m_info = InvalidInput;
		return;
	}

	this->initFactorization(a);

	int info = 0;
	RealScalar recip_pivot_growth, rcond;

	StatInit(&m_sluStat);
	SuperLU_gsisx(&m_sluOptions,
				  &m_sluA,
				  m_q.data(),
				  m_p.data(),
				  &m_sluEtree[0],
				  &m_sluEqued,
				  &m_sluRscale[0],
				  &m_sluCscale[0],
				  &m_sluL,
				  &m_sluU,
				  NULL,
				  0,
				  &m_sluB,
				  &m_sluX,
				  &recip_pivot_growth,
				  &rcond,
				  &m_sluStat,
				  &info,
				  Scalar());
	StatFree(&m_sluStat);

	// FIXME how to better check for errors ???
	m_info = info == 0 ? Success : NumericalIssue;
	m_factorizationIsOk = true;
}

#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename MatrixType>
template<typename Rhs, typename Dest>
void
SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& x) const
{
	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call "
										"either compute() or analyzePattern()/factorize()");

	const int rhsCols = b.cols();
	eigen_assert(m_matrix.rows() == b.rows());

	m_sluOptions.Trans = NOTRANS;
	m_sluOptions.Fact = FACTORED;
	m_sluOptions.IterRefine = NOREFINE;

	m_sluFerr.resize(rhsCols);
	m_sluBerr.resize(rhsCols);

	Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor>> b_ref(b);
	Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor>> x_ref(x);

	m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
	m_sluX = SluMatrix::Map(x_ref.const_cast_derived());

	typename Rhs::PlainObject b_cpy;
	if (m_sluEqued != 'N') {
		b_cpy = b;
		m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
	}

	int info = 0;
	RealScalar recip_pivot_growth, rcond;

	StatInit(&m_sluStat);
	SuperLU_gsisx(&m_sluOptions,
				  &m_sluA,
				  m_q.data(),
				  m_p.data(),
				  &m_sluEtree[0],
				  &m_sluEqued,
				  &m_sluRscale[0],
				  &m_sluCscale[0],
				  &m_sluL,
				  &m_sluU,
				  NULL,
				  0,
				  &m_sluB,
				  &m_sluX,
				  &recip_pivot_growth,
				  &rcond,
				  &m_sluStat,
				  &info,
				  Scalar());
	StatFree(&m_sluStat);

	if (x.derived().data() != x_ref.data())
		x = x_ref;

	m_info = info == 0 ? Success : NumericalIssue;
}
#endif

#endif

} // end namespace Eigen

#endif // EIGEN_SUPERLUSUPPORT_H
